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Abstract 

We consider a mechanism for area preserving Hamiltonian systems which leads to 
the enhanced probabihty, P{X,t), to find small values of the finite time Lyapunov 
exponent, A. In our investigation of chaotic dynamical systems we go beyond the 
linearized stability analysis of nearby divergent trajectories and consider folding of 
the phase space in the course of chaotic evolution. We show that the spectrum of the 
Lyapunov exponents F{X) = limj^oo In -P(A, t) at the origin has a finite value 
F(0) = —A and a slope F'{0) < 1. This means that all negative moments of the 
distribution (e""*'^*) are saturated by rare events with A — s- 0. Extensive numerical 
simulations confirm our findings. 
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An exponential divergency of nearby trajectories in phase space is commonly 
considered as a paradigm of classical chaos [1,2]. The time evolution of the 
distance between two trajectories is determined by the stability matrix 
whose largest eigenvalue grows exponentially like e'^*. For finite time t the value 
of A depends on initial conditions and it is called a finite time Lyapunov expo- 
nent. The probability distribution of finite time Lyapunov exponents P{X,t) 
[2,3,4], especially its behavior at small A, is important for many applications, 
where the measured quantity is sensitive to the existence of trajectories stay- 
ing close for anomalously long time. Examples range from the problems of 
ocean acoustics [5] and branching of 2d electron flow [6] to Loschmidt echo [7] 
and mesoscopic superconductivity [8]. 
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An evolution of small areas in phase space in the linearized approximation is 
described by a combination of area preserving stretching and squeezing (we 
assume two canonical variables x and p). In spite of a broad literature in the 
field the research which goes beyond the linearized stability analysis of dynam- 
ics of close trajectories is sparse. For dissipative systems it was done almost 
twenty years ago by Politi et a/. [9] and by Ott et a/. [10]. In those papers the 
authors considered the consequences of tangency between stable and unstable 
manifolds for chaotic attractors of nonhyperbolic two-dimensional maps. We 
are not aware that similar analysis was applied to the Hamiltonian systems. 
In this letter we try to fill this gap and consider the effect of bending of phase 
space on P(A, t) for area preserving systems. Although our approach is sim- 
ilar to the approach of refs. [9,10], their results were derived for the fractal 
dimension spectrum rather than for the finite time Lyapunov spectrum. Both 
characteristics are related to each other for dissipative systems (see Grass- 
berger et al. in ref.[3]). However for Hamiltonian systems only the Lyapunov 
spectrum is meaningful. Its behavior at small A, which is important for various 
practical applications [5,6,7,8], remains not covered in the existing literature. 

As we will see, a creation of narrow folds in the course of mapping of the phase 
plane results in a great enhancement of the probability to find small values 
of A. As a result all negative moments (e""*^*) for sufficiently large times are 
saturated by the rare events with A ^ 0. 



For large times the probability distribution of finite-time Lyapunov exponents 
P(A, t) [2,3,4] has a generic form (A > 0) 



with a model specific function F{X) called spectrum [3]. The conventional 
selfaveraging Lyapunov exponent Aq appears in the expansion of spectrum 
around the maximum F — (A— Ao)^/2A2 (obviously (A) = Aq, varA = A2/Hor 
Xot ^ 1). If the evolution of the system is described by a Markovian process, 
i.e. if P(A, t) for large t is given by a convolution of P(A, ti) for smaller time 
intervals, then naturally F{0) = —oo. On the other hand, in case of mixed 
phase space one has F(0) = 0, due to the possibility for chaotic orbits to get 
"stuck" on the KAM tori [3,11]. In this letter we consider the distribution 
of Lyapunov exponents in the different regime of fully developed dynamical 
chaos (i.e. in the absence of stable periodic orbits). We argue that in this case 
the spectrum still has a finite limit F{0) = const < 0. 

Dynamics of close trajectories is determined by the stability matrix via 



P(A,i)~exp[iF(A)], 



(1) 




■,M 




) 



(2) 



2 





-0.2 
-0.4 

-0.6 
-0.8 
-1 
-1.2 
-1.4 





1 / 1 1 1 

^ § ■ 


• OOb / rtOB 
• OOm / <>t)B 




' '"o. As." 


i 

■ 


•OO ■ /oo_" 

-"•A:- II 

■ °" 1 1 


■ o 

@ * 

o 

1 1 1 1 



0.5 



1.5 
X 



2.5 



Fig. 1. The finite-time spectrum Ft{\) for the kicked rotator for different kicking 
strengths K = 11(1) and K = 20(11), and times t = 10(b), 15(o), 20(o), 25(»), and 
for the time dependent random potential for two regimes a(v) and b(A) described 
in the text. In case b due to the small values of Lyapunov exponent both i^t(A) and 
A were multiplied by 5. The solid line shows a slope F/(A) = 1. 
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Fig. 2. Distribution of the endpoints of trajectories (xt,pt) having a small value of 
LE A < Ao/18, Aq = ln(iir/2) for the kicked rotator with K = 20 after 5 iterations. 

The finite time Lyapunov exponent, specific for a given trajectory, is defined 
as [12] 

' ' (3) 



A = {2t)-^ In tTM^M/2 + ^J{tiM^M/2f 



Before presenting the analytical theory, we show in Fig. 1 the function Ft{X) = 
ln(P(A, t))/t, found numerically for two models. The first model, which we will 
primarily use, is the kicked rotator [1]: 



Pn+l ^Pn + K sin Xn, X^+l ^ + Pn+l- 



(4) 



The spectrum F{X) in Eq. (1) is strictly defined in the limit t oo. Due to 
A > 0, its finite time approximations Ft{X) always have Ft(0) = —oo, but this 
singularity becomes weaker as time increases. This tendency can be observed 
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Fig. 3. Schematic presentation of the evolution of a small phase space area, described 
in the text. 

in Fig. 1, where different finite time plots are presented (see especially the 
case for K = 11). Data on the figure includes ~ 10^^ trajectories. Calculation 
of the spectrum for longer times would require the unrealistic increase of the 
statistics. However the exponential divergency of close trajectories is already 
very noticeable, being 5x{t)/5x{0) ^ e^^ for largest times in Fig. 1. 

The second model describes a particle moving under the action of a time 
dependent random force, which we introduced via the "random kick" Hamil- 
tonian 

2 N 

H ^^ + KY,Kt-n) sin qx + cos qx) (5) 

2 n q=l 

with random coefficients la^""-*!, \b^q^\ < 1. Both plots in Fig. 1 are for k = 1: 
the plot (a) has = 3 and 50 random kicks, while the plot (6) has N = 1 and 
300 kicks. The behavior of the spectrum for large A is completely different for 
the two models. However, a shape of F{\) at small A appears to be surprisingly 
similar. The results of all simulations are consistent with a linear increase at 
A < Ao, i.e. F(0) = const, F'(0) = 1. 

The distribution of finite time Lyapunov exponents, Eq. (1), is a rough statis- 
tical characteristic which ignores correlations of values of A between different 
trajectories. Such correlations are seen in Fig. 2, where we show the endpoints 
of trajectories, which have a small A. The points on the phase plane are not 
distributed uniformly, as it would be for a random process, but form a pattern 
of narrow curved lines. Such a structure can be explained by the mechanism, 
schematically presented in Fig. 3. The figure shows a time evolution of small 
area of the phase plane around some classical trajectory {a;(r),p(r)}. This evo- 
lution can be conveniently described by the transformation of a small piece 
of square lattice covering the area (Fig. 3a). Few stages of such an evolution 
are shown. (Estimates, analogous to that of Eqs. (6-12) below, were used in 
Refs. [9,10] in order to describe the fractal dimension spectrum /(a) of chaotic 
attr actor.) 

The first stage, shown in Fig. 3b, consists of linear extension and squeezing. 
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After the appropriate orientation of coordinate frame {x, p ^ u,v) it may be 
expressed as 

u Me^»*i , V ve-^°^K (6) 

At this moment we ignore the fluctuations of the Lyapunov exponents and use 
Ao for estimates. The hnear transformation (6) always describes the evolution 
of infinitesimally small areas. For small but finite regions of phase space a 
nonlinear term should be added to Eq. (6) 

u ^ u , V ^ V + (7) 

which is shown in Fig. 3c (the orientation of initial axes (a) is different for 
Figs. 3b, c and d). Parabolas in Fig. 3c have a curvature comparable to the 
size of the system, or the typical momentum, both assumed to be of order 
unity. A consequent evolution of the area consists of further stretching and 
squeezing. However, the map, built from straight lines and parabolas shown 
in Fig. 3c, is robust under this transformation. Squeezing of a parabola gives 
another, more narrow, parabola. It is shown in Fig. 3d and may be expressed 
as 

u ue~^°*^ , V ve^°*^. (8) 

Since the orientation of unstable manifold changes by tt along parabola, the 
squeezing always leads to diminishing of the exponential divergency in some 
region near the vertex of the parabola. Thus a stripe with small values of A is 
formed along the centres of a set of narrow parabolas shown in Fig. 3d. The 
curves in Fig. 2 show these stripes of small A. Combining together eqs. (6-8) 
we express new coordinates u and v through initial ones uq and vq 

u = iioe^°(*i-*^) , V = ^;oe^°(*^-*i) + !ie^°(*^+2*i) . (9) 

Now the Lyapunov exponent is calculated with the help of the stability matrix 
M viaARilnM|/2t (3): 

A = - In [2ch(2Ao(ti - ts)) + ii'e'^o(*2+*)] , (10) 

where t = ti + t2. Thus the finite time Lyapunov exponent depends only on 
one variable u. The function \{u) (10) has a narrow minimum with the depth 
^min — Ao|ti — ^2!/^, and the width 5u = exp(Ao(|ti — — t — ^2))- This 
explains the pattern seen in Fig. 2. Away from the minimum we write 

A = Ao + (ln|M| + Aoi2)A or \u\ = e(^-^»)*-^o*2_ ^^^^ 

We see that at the time t small values of A are found in narrow stripes on the 
plane, each corresponding to the set of a squeezed parabolas (Fig. 3d). 
Among these lines those created at ti ~ t/2 have values as small as A ~ 1/t, 
leading to nonvanishing P(A, t) at the origin: F(0) > —00. In order to estimate 
the probability distribution P{\,t) at A <^ Ao, we first notice that a fold 
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Fig. 4. nX versus n for kicked rotator witli K = 20. Examples of typical trajectories 
(slope ln(iC/2), shown dashed), and rare trajectories having A < 1/5 after 15 
iterations (shown by solid lines). 




Fig. 5. The contour plot of the joint distribution W{\ Ini?) of the Lyapunov expo- 
nent and the logarithm of curvature, Eq. (19), for the kicked rotator with K = 
after 15 iterations. Contours were drawn by taking 10^^ initial conditions and al- 
locating the results into the bins of two-dimensional histogram. The number of 
counts increases by 10 for each black contour, from outer to inner (the probability 
on the inner border of the contour is twice the probability on the outer side). The 
correlation between small A and small R confirms our phase-space-folding induced 
mechanism. 

created at time ti had a length ~ 1. It becomes as long as e'*'''*^ at the time t. 
Since the probability to find the Lyapunov exponent in the interval d\ is given 
by the area, where such Lyapunov exponents are found, we may write 

P(A, t) ~ e^'^^'du/dX ~ exp[(A - Xo)t]. (12) 

We keep only exponentially large contributions in this estimate. 

In derivation of (12) we ignored fluctuations of Lyapunov exponents during 
the squeezing stages (6) and (8). Due to fluctuations P{X,t) is saturated by 
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only exponentially small fraction of strongest folds. Still wc may say that the 
simple folding mechanism described by the eqs. (6-9) is sufficient for having 
a finite value of the spectrum at the origin F(0) > — oo. It is also capable to 
prevent the probabihty from decaying at small A faster than P{X,t) cx c^*. 
Thus we write 

^(0) = -A > -Ao , F'(0) < 1. (13) 

Taking into account fluctuations of the Lyapnov exponent during the stages 
(6) and (8) may be done as 

P{X,t)^ I e^*-^i*i-^^*^P(Ai,ii)P(A2,i2) 

X 6{ti +t2- t)6{Xiti - X2t2 - Xt)dtidt2dXidX2 . (14) 

Detailed investigation of this formula goes beyond the scope of this paper. The 
behavior P(A, t) at small A may however be analyzed easily. One readily shows, 
that the spectrum of the form (13) with F' < is consistent with the eq. (14). 
Substitution of P{X,t) in the form (13) into the r.h.s. of (14) leads to the 
same small A behavior in the l.h.s. This docs not happen for F'{0) > 1. Also 
if F' < 0, the integral over Ai and A2 for small A is saturated by the region 
Ai,A2 <^ Aq. Figure 4 shows nX ~ \n{Sxn/SxQ) for few typical trajectories 
and for the rare trajectories with anomalously small value of A. As we see, 
typically two initially very close trajectories would escape from each other 
almost without fluctuations 6x(t) ~ e^"^dx{0). For rare trajectories, leading 
after large time t to X(t) ~ 0, the extent of linear divergency is always much 
smaller, 5x{t') -C e'^"^' 5x{0) for any t' < t, in agreement with (14). 

Although our result eq. (13) gives only a bound on how fast the spectrum may 
increase with A , F'{0) < 1, numerical simulations of Fig. 1 are consistent with 
F'(0) ~ 1. Even in this "weak" form the result (13) has a strong predictive 
power, leading to (m = 1, 2, ...) 

lim ln(e-™^*) = -A. (15) 

t— »oo 

The dominant contribution to the averaged value here comes from the rare 
trajectories with e~^* ~ 1. An exponentially small value of (e~"^'^*) is due to 
a small number of such trajectories. Eq. (15) should be compared with the 
often used Gaussian approximation ln(e~"^'^*)G ~ —m{Xo — 111X2/2). 

Folding of the phase plane requires an introduction of a new quantitative 
characteristic, in addition to Lyapunov exponent A describing the linearized 

divergency of the trajectories. To describe folding we introduce the radius of 
curvature i? of a small area bended in the course of evolution. To this end, we 
first analyze a nonlinear divergency of the closed trajectories by expanding 
the map (4) to the second order in the displacement between trajectories. If 
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we write two components of the displacement 5p, 5x as a column q, then 



qo Aqo + M ^ I qo ^nqo, 



(16) 



where M.n-,^n)^n are 2x2 matrices. The stabiUty matrix for kicked rotator 
is given by a product 



„ , 1 Kcosxj 

Mn=\{ 

i<n \ 1 1 + K COS Xi 



(17) 



Two other symmetric matrices may be found from the recursion relations [13]. 



An+l = Ai + (1 + K COS Xn)Bn - KsiuXnCn, 

Bn+1 ^ An+1 + Bn, (18) 

where Ai — Bi — and the elements of the matrix C„ are {Cn)ij = {M.n)i2{M.n)j2- 
Consider the radius of curvature R of the image of a line = (0, Sx) as a 
measure of bending (similar results may be obtained for the curvature of the 
image of line qo = {5p, 0)). A simple calculation gives 

1 ^ {An)22{Mn)22 - {Bn)22{Mn)l2 . ^ 

In order to analyze correlations between folding and the values of the finite 
time Lyapunov exponent we introduce the joint distribution 1^(A, Ini?). It is 
shown in Fig. 5 [14] . We see that for given A the curvature may vary by orders 
of magnitude. However the general trend is transparent: the smaller values of 
A correspond to the smaller R [15]. 

The finite fraction of trajectories with A ~ in our mechanism does not 
require the existence of any stable islands. Still our main result (13) includes 
also the case of mixed phase space, corresponding to A = 0. In the presence of 
very small stable islands one expects a competition of two effects, when our 
mechanism will govern the small A spectrum at intermediate times, while the 
limit t — > oo will be due to dynamical traps [11] in the mixed phase space. 

In conclusion, for the case of developed chaos, the physical quantities which 
may be calculated or measured naturally have a statistical meaning. In this 
letter we have shown, how in such systems folding of the x, p-plane inevitably 
leads to creation of elongated areas, where initially diverging trajectories re- 
vert to approaching to each other again. This allowed us to find a uniform 
description of the tail of distribution of finite time Lyapunov exponents, and 
to calculate all the negative moments of the distribution (15). Negative here 
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means that wc are able to describe the physical quantities which are sensitive 
to the contributions where, in-spite of chaos, trajectories stay longer together 
and do not diverge. 

One immediate application of our findings will be a calculation of the Loschmidt 
echo [7], which was shown [16] to decay on average like {e~^^). This behavior 
should now be reconsidered in view of our eq. (15). Another problem of cur- 
rent interest, where the result is sensitive to the existence of long trajectories 
with small A, is a calculation of the gap in the spectrum of the Andreev bil- 
liard [8]. Application of our result Eq. (13) to this problem leads to the strong 
modification of the low energy Andreev spectrum and substantial lowering of 
the superconducting gap [17]. 

Discussions with C. W. J. Beenakker, A. N. Morozov, J. Tworzydlo, and 
J. Whitakcr are greatly appreciated. This work was supported by the Dutch 
Science Foundation NWO/FOM and by the SFB TR 12. 
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